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ABSTRACT 

Strong gravitational lenses with measured time delays between the multiple images and models of 
the lens mass distribution allow a one-step determination of the time-delay distance, and thus a mea- 
sure of cosmological parameters. We present a blind analysis of the gravitational lens RXJ1131 — 1231 
incorporating (1) the newly measured time delays from COSMOGRAIL, the COSmological MOnitor- 
ing of GRAvItational Lenses, (2) archival Hubble Space Telescope imaging of the lens system, (3) a new 
velocity-dispersion measurement of the lens galaxy of 323 ±20 kms^^ based on Keck spectroscopy, and 
(4) a characterization of the linc-of-sight structures via observations of the lens' environment and ray 
tracing through the Millennium Simulation. Our blind analysis is designed to prevent experimenter 
bias. The joint analysis of the data sets allows a time-delay distance measurement to 6% precision that 
takes into account all known systematic uncertainties. In combination with the Wilkinson Microwave 
Anisotropy Probe seven-year (WMAP7) data set in flat wCDM cosmology, our unblindcd cosmological 
constraints for RXJ1131-1231 are: Hq = 80.0l^;^kms-i Mpc-\ n^e = 0.79± 0.03, w = -1.25to ll- 
We find the results to be statistically consistent with those from the analysis of the gravitational lens 
B1608+656, permitting us to combine the inferences from these two lenses. The joint constraints from 
the two lenses and WMAP7 arc Hq = 75.2+4;2kms-i Mpc-\ Qdc = 0.76+o:o3 and w = -1.14^o;2o in 
flat wCDM, and Hq = 73.1^3:6 km s^^ Mpc~\ = 0.75tg;^2 and f^k = 0.003lJ^;^^^ in open ACDM. 
Time-delay lenses constrain especially tightly the Hubble constant Hq (5.7% and 4.0% respectively in 
wCDM and open ACDM) and curvature of the universe. The overall information content is similar 
to that of Baryon Acoustic Oscillation experiments. Thus, they complement well other cosmological 
probes, and provide an independent check of unknown systematics. Our measurement of the Hubble 
constant is completely independent of those based on the local distance ladder method, providing an 
important consistency check of the standard cosmological model and of general relativity. 
Subject headings: galaxies: individual (RXJ1131 — 1231) — gravitational lensing: strong — methods: 
data analysis — distance scale 



L INTRODUCTION 

In the past century precise astrophysical measurements 
of the geometry and content of the universe (hereafter 
cosmography) have led to some of the most remarkable 
discoveries in all of physics. These include the expansion 
and acceleration of the Universe, its large scale structure, 
and the existence of non-baryonic dark matter (see re- 
view bv lFreedman fc Turned f200l . These observations 
form the empirical foundations of the standard cosmo- 
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logical model, which is based on general relativity and 
the standard model of particle physics but requires addi- 
tional non-standard features such as non-baryonic dark 
matter and dark energy. 

Even in the present era of so-called precision cosmogra- 
phy, many profound questions about the Universe remain 
unanswered. What is the nature of dark energy? What 
are the properties of the dark matter particle? How many 
families of rclativistic particles are there? What are the 
masses of the neutrinos? Is general relativity the correct 
theory of gravity? Did the Universe undergo an infla- 
tionary phase in its early stages? 

From an empirical point of view, the way to address 
these questions is to increase the accuracy and precision 
of cosmographic experiments. For example, clues about 
the nature of dark energy can be gathered by measur- 
ing the expansion history of the Universe to very high 
precision, and modeling the expansion as being due to 
a dark energy component having an equation of state 
parametrized by w that evolves with cosmic time (e.g., 
Fricman ct al. 200^, and references therein). Likewise, 
competing inflationary models can be tested by measur- 
ing the curvature of the Universe to very high precision. 
Given the high stakes involved, it is essential to develop 
multiple independent methods as a way to control for 
known systematic uncertainties, uncover new ones, and 
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ultimately discover discrepancies that may reveal new 
fundamental physics. For example, a proven inconsis- 
tency between inferences at high redshift from the study 
of the cosmic microwave background, with inferences at 
lower redshift from galaxy redshift surveys would chal- 
lenge the standard description of the evolution of the 
Universe over this redshift interval, and possibly lead to 
revisions of either our theory of gravity or of our assump- 
tions about the nature of dark matter and dark energy. 

In this paper we present new results from an obser- 
vational program aimed at precision cosmography using 
gravitational lens time delays. The idea of doing cos- 
mography with time-de lay lenses goe s back fifty years 
and it is a simple one (IRefsdall[T96l . When a source 
is observed through a strong gravitational lens, multiple 
images form at the extrema of the time-delay surface, ac- 
cording to Fcr mat's principle (e.g. . [Schneider et al.lll992l : 
iFalcQ 2005: Schneider et al.ll2006f ). If the source is vari- 
able, the time delays between the images can be mea- 
sured by c areful rnqnitor ing of the image light curves 
(see, e.g., iCourbinI I2003D . With an accurate model 
of the gravitational lens, the absolute time delays can 
be used to convert angles on the sky into an abso- 
lute distance, the so-called time-delay distance, which 
can be compared with predictions from the cosmolog- 
ical model given the lens and source re d shifts (e.g. , 
IBlandford fc NaravanI [1991 iJacksonI [20071: [Si3 iMfl 
and references therein). This distance is a combination 
of three angular diameter distances, and so is primar- 
ily sensitive to the Hubble constant, with some higher 
order dependence on the other cosmolo gical parameters 
(|Coe fc MoustakasI l2009t iLindeii |2011[ ). Gravitational 
time delays are a one-step cosmological method to de- 
termine the Hubble constant that is completely in depen- 
dent of the local cosmic distance ladder ( ^ecdman et al 
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iRiess et al.ll2011l : lFreedman et al.l[2012l : iReid et"aT 



20121) . Knowledge of the Hubble constant is currently 
the key limiting factor in measuring parameters like the 
dark energy equation of state, curvature, or neutrino 
mass, in combination with other probes like the cos- 
mic microwave background (iFre cdman fc Madord 120101: 
Riess et all 120111: iFreedman et al...2012; iWeinberg et al.l 



2012t iSuvu et all I2012D . These features make strong 



gravitational time delays a very attractive probe of cos- 
mology. 

Like most high-precision measurements, however, a 
good idea is only the starting point. A substantial 
amount of effort and observational resources needs to be 
invested to control the systematic errors. In the case of 
gravitational time delays, this has required several ob- 
servational and modeling breakthroughs. Accurate, long 
duration, and well-sampled light curves are necessary to 
obtain accurate time delays in the presence of microlcns- 
ing. Modern light curves have mu ch higher photometric 
precision, sanipling and duration (jFassnacht et al.|[2002l : 
iCourbin et al.l l201ll) compared to the early pioneering 
light curves (e.g.. lLehar et al.|[l992| ). High resolution im- 
ages of extended features in the source, and stellar kine- 
matics of the main deflector, provide hundreds to thou- 
sands of data points to constrain the mass model of the 
main deflector, thus reducing the degeneracy between the 
distance and the gravitational potential of the lens that 
affected previous models constrained on ly by the pos i- 
tions of the lensed quasars (e.g., .Schechter et al.|[l997D . 



Finally, cosmological numerical simulations can now be 
used to characteriz e the distribution o f mass along the 
line of sight (LOS) (iHilbert et al.ll2009f ). which was usu- 
ally neglected in early studies that were not aiming for 
precisions of a few percent. The advances in the use of 
gravitational time delays as a cosmographic probe are 
summarized in the a nalysis of the grav itational lens sys- 
tem B1608-h656 by iSuvu et all (|2010[) . In that paper, 
we demonstrated that, with sufficient ancillary data, a 
single gravitational lens can yield a time-delay distance 
measured to 5% precision, and the Hubble constant to 
7% precision. In combination with the Wilkinson Mi- 
crowave Anisotropy Probe 5-year (WMAP5) results, the 
B1608-I-656 time-delay distance constrained w to 18% 
precision and the curvature parameter to ±0.02 preci- 
sion, comparable to c ontemporary Baryon Acoustic Os- 
cillation experiments (jPercival et al.l[2"007j ) and observa- 
tions of the growth of mass ive galaxy clusters (on w con- 
straints; iMantz et ai]|2010[) . 

Building on these recent developments in the analysis, 
and on the state-of-the-art monitoring campaigns carried 
out by the COSMOGRAIL (COSmological MOnitoring 
of GRAvItati on al Le nses : e.g.. [Vuissoz et al, ,2QQS 
Courbin et all [20Tlt iTewes et all l2012bh and 
Kochanek et al.l (l2006f ) teams, it is now possible to 



take gravitational time delay lens cosmography to the 
next level and achieve precision comparable to current 
measurements of the Hubble cons tant, flatness, w 
and other cosm ological parameters ([Riess et al.l 120111 : 
IFreedman et al.ii2012; Komatsu et al.ll2011[ ). To this end 
we have recently initiated a program to obtain data and 
model four additional gravitational lens systems with 
the same quality as that of B1608+656. 

We present here the results for the first of these sys- 
tems, RXJ1131— 1231, based on new time delays mea- 
sured by the COSMOGRAIL coUaboration (jTewes et al.l 
l2012bD . new spectroscopic data from the Keck Tele- 
scope, a new analysis of archival Hubble Space Telescope 
(HST) images, and a characterization of the LOS ef- 
fects through numerical simulations, the observed galaxy 
number counts in the field, and the modeled external 
shear. We carry out a self-consistent modeling of all 
the available data sets in a Bayesian framework, and in- 
fer (1) a likelihood function for the time-delay distance 
that can be combined with any other independent probe 
of cosmology, and (2) in combination with our previ- 
ous measurement of B1608+656 and the WMAP 7-year 
(WMAP7) results, the posterior probability density func- 
tion (PDF) for the Hubble constant, curvature density 
parameter and dark energy equation-of-state parameter 
w. 

Three additional lens systems are scheduled to be ob- 
served with HST in cycle 20 (GO 12889; PI Suyu) and 
will be published in forthcoming papers. An integral part 
of this program is the use of blind analysis, to uncover 
unknown systematic errors and to avoid unconscious ex- 
perimenter bias. Only when each system's analysis has 
been judged to be complete and final by its authors, are 
the implications for cosmology revealed. These results 
are then published without any further modification. In 
this way, we can assess whether the results are mutually 
consistent within the estimated errors or whether un- 
known systematics are adding significantly to the total 
error budget. 
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This paper is organized as follows. After a brief re- 
cap of the theory behind time-delay lens cosmography in 
Section[2l we summarize our strategy in Section[3]and de- 
scribe our observational data in Section HI In Section [5l 
we write out the probability theory used in the data 
modeling and describe the procedure for carrying out the 
blind analysis. The lensing and time-delay analysis are 
presented in Section|6l and a description of our treatment 
of the LOS mass structure in the RXJ1131— 1231 field 
is in Section [71 We present measurements of the time- 
delay distance, and discuss the sources of uncertainties in 
Section [HI We show our unblinded cosmological parame- 
ter inferences in Section [HI which includes joint analysis 
with our previous lens data set and with WMAP7. Fi- 
nally, we conclude in Section [TOl Throughout this paper, 
each quoted parameter estimate is the median of the ap- 
propriate one-dimensional marginalized posterior PDF, 
with the quoted uncertainties showing, unless otherwise 
stated, the 16**^ and 84*'^ percentiles (that is, the bounds 
of a 68% credible interval). 

2. COSMOGRAPHY FROM GRAVITATIONAL LENS TIME 
DELAYS 

In this section, we give a brief overview of the use 
of gravitational lens time delays to study cosmology. 
More details on the subject c an b e fouiid in, e.g., 
Schneider etldl ([20061) . [Jacksonl (|2007f ) . (fMO ) and 
Suvu et al.l ( 2010( r Readers familiar with time-delay 
lenses may wish to proceed directly to Section [3] 

In a gravitational lens system, the time it takes the 
light from the source to reach us depends on both the 
path of the light ray and also the gravitational potential 
of the lens. The excess time delay of an image at angular 
position 9 = {9 1,62) with corresponding source position 
(3 = ^2) relative to the case of no lensing is 



t(0,/3) 



Dm 



{6 -(if 



(1) 



where Da* is the so-called time-delay distance, c is the 
speed of light, and is the lens potential. The time- 
delay distance is a combination of the angular diameter 
distance to the lens (or deflector) {D^) at redshift Zd, 
to the source (-Dg), and between the lens and the source 

DAt^il + z^)^. (2) 
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The lens potential ip{9) is related to the dimensionless 
surface mass density of the lens, k{9), via 



V^xjj{9) = 2k{9), 



where 



49) 



^{Dd9) 



-'crit 
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Y,{Dd9) is the surface mass density of the lens (the pro- 
jection of the three dimensional density p along the LOS), 
Ecrit is the critical surface mass density defined by 



Ecrit 



and G is the gravitational constant. 



(5) 



For lens systems whose sources vary in time (as do ac- 
tive galactic nuclei, AGNs), one can monitor the bright- 
nesses of the lensed images over time and hence measure 
the time delay, Aty , between the images at positions 9i 
and 9i: 



AU,=t(9,,f3)-t{9j,f3) 



DAt 



X0,)(6) 



By using the image configuration and morphology, one 
can model the mass distribution of the lens to determine 
the lens potential ^{9) and the unlensed source position 
/3. Lens systems with time delays can therefore be used 
to measure D^t via Equation ([6]), and constrain cosmo- 
logi cal models via the distance-redshift tes t (e.g . jRefsdail 
[1^4, 196a iFadelv ern][2nTnl : [Suvu et al.l[2nT(Tf ). Having 
dimensions of distance. Da* is inversely proportional to 
iJo, and being a combination of three angular diameter 
distances, it depends weakly on the other cosmological 
parameters as well. 

The radial slope of the lens mass distribution and the 
time-delay distance both have direct influence on the ob- 
servables: for a given time delay, a galaxy with a steep 
radial profile leads to a lowe r Dai than that of a galaxy 
with a shallow profile (e.g., iWitt et al.ll2000l: iWucknitd 
120021: lKochanekl[2002| ). Therefore, to measure Da*, it 
is necessary to determine the radial slope of the lens 
galaxy. Several authors have shown that the spatially 
extended sources (such as the host galaxy of the AGN 
in time-delay lenses) can be used to measure the ra- 
dial slope a t the image positions, wher e it matters (e.g. , 
Dve fc Warr en 2005: Dve et al. 200^ iSuvu et al.l I2OIOI : 



Vegetti et al.]l20ldHSuvTj[2012t) . 



In addition to the mass distribution associated with 
the lens galaxy, structures along the LOS also affect the 
observed time delays. The external masses and voids 
cause additional focussing and defocussing of the light 
rays respectively, and therefo re affect the time delays and 
Dm. i nferences. We follow [Keetonl ()2003l ). iSuvu et aD 
(|2010( ) and many others and suppose that the effect of 
the LOS structures can be characterized by a single pa- 
rameter, the external convergence Atext, with positive val- 
ues associated with overdense LOS and negative values 
with underdense LOS. Except for galaxies very nearby 
to the strong lens system, the Kcxt contribution of the 
LOS structures to the lens is effectively constant across 
the scale of the lens system. 

Given the measured delays between the images of a 
strong lens, a mass model that does not account for the 
external convergence leads to an under/overprediction of 
Dai for over/underdense LOS. In particular, the true 
Dai is related to the modeled one by 



Dai 



remodel 
1 ^ ^nxt 



(7) 



Two practical approaches to overcome this degener- 
acy are (1) to use the stel lar kinemat ics of the lens 
galax y (e.g. .iTreu fc Koopma ns 2002: Ko opmans fc Treul 
200I [iVeu fc Koopman'g~l2004; .Barnabe et al.l 120091: 



Auger et al l 120101 : ISuvu et al. \ \20m ISonnenfeld et al.l 
201a) " to make an independent estimate of the lens mass 
and (2) to study the environment of the lens system 
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e.g.. iKeeton fc Zabludofl [20p: iFassnacht et all [20061: 



Momcheva et all 120061: iSuvu et al .l l2010t iW ong et al l 
201lHFassnacht et al.ll20f¥ in order to estimate K^xt 

di- 
rectly. In Section [21 we combine both approaches to infer 

3. ACCURATE AND PRECISE DISTANCE 
MEASUREMENTS 

We summarize our strategy for accurate and precise 
cosmography with all known sources of systematic un- 
certainty taken into account. We assemble the following 
key ingredients for obtaining Z^At via Equation ^ 

• observed time delays: dedicated and long-duration 
monitoring, particularly from COSMOGRAIL, 
yield s delays with uncer tainties of only a few per- 
cent (|Tewes et al.ll2012al [bl). 

• lens mass model: deep and high-resolution imag- 
ings of the lensed arcs, together with our flexible 
modeling techniques that use as data the thousands 
of surface brightness pixels of the lensed source, al- 
low constraints of the potential difference between 
the lensed imag es (in Equation (pi )) at the few per- 
cent level (e.g.. ISuvu et al.ll2010r 

• external convergence: the stellar velocity disper- 
sion of the lens galaxy provides constraints on both 
the lens mass distribution and external conver- 
gence. We further calibrate observations of galaxy 
counts in the fields of lenses (jFassnacht et al.ll2011[) 
with ray tracing through n umerical simulation s of 
large-scale structure (e.g.. iHilbert et all I2007D to 
const rain directly and s tatistically /text at the ~5% 
level (|Suvu et aU\2(M ). 

With all these data sets for the time-delay lenses, we 
can measure D^t for each lens with ~ 5 — 8% precision 
(including all sources of known uncertainty). A com- 
parison of a sample of lenses will allow us to test for 
residual systematic effects, if they are present. With sys- 
tematics under control, we can combine the individual 
distance measurements to infer global properties of cos- 
mology since the gravitational lenses are independent of 
one another. 

4. OBSERVATIONS OF RXJ1131-1231 

The gravitat i onal l ens RXJ1131 — 1231 was discovered 
by iSluse et all ()2003() during polarimetric imaging of a 
sample of radio quasars. The spectroscopic redshifts of 
the lens and the quasar source are Zd — 0.295 and Zg = 
0.658, respectively. We present the archival HST images 
in Section [4Jl the time delays from COSMOGRAIL in 
Section 14.21 the lens velocity dispersion in Section 14.3! 
and information on the lens environment in Section 14.41 

4.1. Archival HST imaging 

HST Advanced Camera for Surveys (ACS) images were 
obtained for RXJ1131-1231 in two filters, F814W and 
F555W (ID 9744; PI Kochanek). In each fiher, five expo- 
sures were taken with a total exposure time of 1980s. We 
show in Figure m the F814W image of the lens system. 
The background quasar source is lensed into four images 
denoted by A, B, C, and D, and the spectacular features 
surrounding the quasar images are the lensed images of 



B 




D 



N 



Fig. 1.— HST ACS F814W image of the gravitational lens 
RXJ1131 — 1231. The lensed AGN images of the spiral source 
galaxy are marked by A, B, C and D, and the star forming re- 
gions of the spiral galaxy form the spectacular lensed structures. 
The primary lens galaxy and the satellite lens galaxy are indicated 
by G and S, respectively. 



the q uasar host that is a spiral galaxy (jClaeskens et al.l 
l2Q06f ). The primary lens galaxy is marked by G, and 
the object marked by S is most likely a satellite of G 
(jClaeskens et al.ll2006[ ). Henceforth, we refer to S as the 
satellite. 

We reduce the images using MultiDrizzle[3 with 
charge transfer inefficiency taken into acc ount (e.g., 
[Anderson fc BedinI [MotlMassev et~alll2010[) . The im- 
ages are drizzled to a final pixel scale of 0'.'05 pix~^ and 
the uncertainty on the flux in each pixel is estimated from 
the science and the weight image by adding in quadra- 
ture the Poisson noise from the source and the back- 
ground noise due to the sky and detector readout. We 
note that in some of the exposures the central regions of 
the two brightest AGN images arc slightly saturated and 
are masked during the drizzling process. 

To model the lens system using the spatially extended 
Einstein ring of the host galaxy, we focus on the F814W 
image since the contrast between the ring and the AGN 
is more favorable in F814W. In particular, the bright 
AGNs in F555W have diffraction spikes extending into 
the Einstein ring that are difficult to model and are thus 
prone to systematics effects. The detailed modeling of 
the F814W image is in Section [HI 

4.2. Time delays 

We use the new tim e-delay mea s uremen ts of 
RXJ1131-1231 presented in iTewes et all (1201 2bf ). The 
COSMOGRAIL and Kochanek et al. teams have moni- 
tored RXJ1131 — 1231 since December 2003 using several 
optical 1-1. 5m telescopes. Resolved light curves of the 
four AGN images are extracted from thes e observations 
by "deconvolution photometry" , following IMagain et al.l 

MultiDrizzle is a product of the Space Telescope Science 
Institute, which is operated by AURA for NASA. 
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()1998D . These curves presently span 9 years with over 
700 epochs, and display a typical sampling of 2-3 days 
within the observation seasons. The time delays are mea- 
sured throu gh several new and independent techniques 
(detailed in iTewes et alJ l2012al ) , all specifically devel- 
oped to handle microlensing variability due to stars in 
the lens galaxy. All these techniques yield consistent re- 
sults, attributed to both the long light curves and the 
comprehensive uncertainty estimation. For our analy- 
sis we select the time-delay measurements from the re- 
gression difference technique, which displays the smallest 
systematic error when applied to synthetic curves mim- 
icking the microlensing variability in RXJ1131— 1231. In 
particular we use the time delays relative to image B, 
namely: AtAB = 0.7 ±1.4 days, AtcB = -0.4 ±2.0 days, 
and A^DB = 91.4 ± 1.5 days, where the uncertainties are 
conservative and direct sums of th e estimated s t atistica l 
and systematic contributions from lTewes et al.1 ()2012bl ). 

4.3. Lens velocity dispersion 

We observed RXJ1131— 12 31 with the Low- Resolution 
Imaging Spectrometer fLRIS: [Oke et al.lll995[ ) on Keck 1 
on 4-5 January 2011. The data were obtained from the 
red side of the spectrograph using the 600/7500 grat- 
ing with the D500 dichroic in place. A slit mask was 
employed to obtain simultaneously spectra for galaxies 
near the lens system. The night was clear with a nomi- 
nal seeing of 0'.'7, and we use 4 exposures of 1200s for a 
total expos ure time of 4800s. 

We follow I Auger et al.l (|2008[ ) to reduce each exposure 
by performing a single resampling of the spectra onto a 
constant wavelength grid. We use the same wavelength 
grid for all exposures to avoid resampling the spectra 
when combining them. An output pixel scale of 0.8 A 
pix"-'^ was used to match the dispersion of the 600/7500 
grating. Individual spectra are extracted from an aper- 
ture 0'.'81 wide (corresponding to 4 pixels on the LRIS red 
side) centered on the lens galaxy. The size of the aperture 
was chosen to avoid contamination from the spectrum of 
the Icnscd AGNs. We combine the extracted spectra 
by clipping the extreme points at each wavelength and 
taking the variance-weighted sum of the remaining data 
points. We repeat the same extraction and coaddition 
scheme for a sky aperture to determine the resolution of 
the output co-added spectrum: R = 2300, corresponding 
to (Tobs = 56 km s~^. The typical signal-to-noise ratio per 
pixel of the final spectrum is ^20. 

The stellar vel ocity dispersion i s determined in the 
same manner as ISuvu et all (l2010l) . Briefly, we use a 
suite of stellar templates of K and G giants, augmented 
with on e A and one F star template, from the INDO-US 
library (jValdes et al.ll200"^ to fit directly to the observed 
spectrum, after convolving each template with a kernel to 
bring them to the same spectral resolution as the data. 
First and second velocity moments are proposed by a 
Markov chain Monte Carlo (MCMC) simulation and the 
templates are shifted and broadened to these moments. 
We then fit the model templates to the data in a lin- 
ear least squares sense, including a fifth order polyno- 
mial to acc ount for any emis sion from the background 
source fe.g. JSuvu et al.ir2010D . The observed and mod- 
eled spectra are shown in Figured! Our estimate for the 
central linc-of-sight velocity dispersion from this infer- 
ence is tr = 323 ± 20kms~^, including systematics from 
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Fig. 2.— Top panel: The LRIS spectrum of RXJ1131-1231 
(black line) with a model generated from 9 INDO-US templates 
and a 5^^ order continuum overplotted (red line, with green show- 
ing the continuum). The gray shaded areas were not included in 
the fit. Bottom panel: the residuals of the model fit. From the 
spectrum and model, we measure a central velocity dispersion of 
a = 323 ± 20kms^^, including systematic uncertainties. 

changing the polynomial order and choosing different fit- 
ting regions. 

4.4. Galaxy counts in the field 

iFassnacht et al.l ()2011|) counted the number of galax- 
ies with F814W magnitudes between 18.5 and 24.5 that 
lie within 45" from the lens system. Compared to the 
aperture counts in random lines of sight in pure-parallel 
fields, RXJ1131— 1231 has 1.4 times the average number 
of galaxy counts. We use this relative galaxy count in 
Section [7^ to estimate statistically the external conver- 
gence. 

5. PROBABILITY THEORY FOR COMBINING MULTIPLE 
DATA SETS 

We now present the mathematical framework for the 
inference of cosmological parameters from the combina- 
tion of the data sets described in the previous section. 
In order to test for the presence of any unknown sys- 
tematic uncertainty, we describe a procedure for blind- 
ing the results during the analysis phase in Section 15.21 
This procedure is designed to ensure against unconscious 
experimenter bias towards "acceptable" results. 

5.1. Joint analysis 

The an alysis performed h ere is similar to the one pre- 
sented in iSuvu et al.l ()20100 . with a few improvements. 
Wc briefly describe the procedure below. 

The data sets are denoted by dACS for the ACS im- 
age (packaged into a vector of 160^ surface brightness 
values), A.t for the delays between the images, a for 
the lens velocity dispersion, and denv for properties of 
the lens environment such as the relative galaxy count 
nr = '^gai/('^gai)- We are interested in obtaining the pos- 
terior PDF of the model parameters ^ given all available 
data, 

-P(€|dACS, At,CT,de„v) OC P(dACS, At,a, de„v|l)P(l) 

where the proportionality follows from Bayes' Theorem. 
The first term is known as the likelihood, and the second 
is the prior PDF. Since the data sets are independent. 
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the likelihood is separable, 

P{d ACS, At , (7 , denv 

P(^|l)P(de„v|0- (9) 

Some of the parameters influence all the predicted 
data sets, while other parameters affect the fitting 
of particular data sets only. Specifically, £ = 
{7r,7',6'E,7ext,''/,?'ani,«;oxt}, where tt are the cosmologi- 
cal parameters (e.g., Hq, w, fide), l' is the radial profile 
slope of the main lens galaxy (where p (x r^^ ), 9-^ is the 
Einstein radius of the main lens (that characterizes the 
normalization of the lens mass profile) , 7oxt is the exter- 
nal shear strength at the lens, r/ denotes the remaining 
lens model parameters for the ACS dataQ ''ani is the 
anisotropy radius for the stellar orbits of the lens galaxy, 
and Koxt is the external convergence. For the Icnsing and 
time delays, we subsume the cosmological dependence 
into the time-delay distance Da* = -DAt(7r). Keeping 
only the direct dependencies in each of the likelihoods, 
we obtain 

-P(l|<iACS, Ai, (7, denv) OC P(<iACS, At|i:)At,7':^E,7ext, 
T/, /Cext)-P(0'|7r, 7', ^E, Tani, ^oxt )^'(denv | ^oxt , 7oxt) 

P(7r)P(7')P(0E)P(7ext)P(j7)-P('-ani)-P(Atext). (10) 

For cosmography, we are interested in the cosmological 
parameters tt after marginalizing over all other parame- 
ters 

P(7r|dACS, At, CT, denv) = / d7' d^E d7oxt d?7 d^ani 

P(||dACS,A*,CT,de„v) (11) 

We describe the forms of the partially marginalized lens- 
ing and time-delay likelihood in Section [51 the kinemat- 
ics likelihood in Section mi and the external convergence 
likelihood in Section 17.21 For marginalizing the param- 
eters that are common to t he data sets, we i mpor tance 
sample the priors following iLewis fc Bridld ()2002l ) and 
iSuvu et all ()2010D (a procedure sometimes referred to as 
"simple Monte Carlo"). 

5.2. Blind analysis 

We blind the analysis to avoid experimenter bias, al- 
lowing us to test for the presence of residual systemat- 
ics in our analysis technique by comparing the final un- 
blinded results from RXJ1131— 1231 with the constraints 
fro m the previous ana lysis of B16084-656. As described 
bv lConlev et al.l (j2006f ). the blinding is not meant to hide 
all information from the experimenter; rather, we blind 
only the parameters that concern the cosmological infer- 
ence. 

We define two analysis phases. During the initial 
"blind" phase, we compute likelihoods and priors, and 
sample the posterior PDF, as given above, taking care 
to only make parameter-space plots using one plotting 
code. This piece of software adds offsets to the cos- 
mological parameters (Da* and the components of tt) 
before displaying the PDFs, such that we always see the 
marginalized distributions with centroids at exactly zero. 
We can therefore still see, and measure, the precision of 

excluding the source surface brightness parameters s that can 
be marginaUzed analytically 



the blinded parameters, and visualize the correlations be- 
tween these parameters, but without being able to see if 
we have "the right answer." Both the parameter uncer- 
tainties and degeneracies serve as useful checks during 
this blind phase: the plotting routine can overlay the 
constraints from different models to investigate sources 
of statistical and systematic uncertainties. 

During the blind phase we performed a number of tests 
on the modeling to quantify the sources of uncertain- 
ties, and to check the robustness of the results. These 
are described in Sections HHSl At the end of the tests, 
the collaboration convened a telecon to unblind the re- 
sults. The authors SHS, MWA, SH, PJM, MT, TT, CDF, 
LVEK, DS, and FC discussed in detail the analysis and 
the blinded results, over a summary website. After all 
agreeing that the blind analysis was complete, and that 
we would publish without modification the results once 
unblinded, a script was run to update automatically the 
same website with plots and tables containing cosmolog- 
ical constraints no longer offset to zero. These are the 
results presented in Sections [8] and 19.11 

6. LENS MODELING 

In this section, we simultaneously model the ACS im- 
ages and the time delays to measure the lens model pa- 
rameters, particularly P'At, 7', ^E, and 7cxt- 

6.1. A comprehensive mass and light model 

The ACS image in Figure [1] shows the light from the 
source as lensed by the galaxies G and S. To predict the 
surface brightness of the pixels on the image, we need 
a model for the lens mass distribution (that deflects the 
light of the source) , the lens light distribution, the source 
light distribution and the point spread function (PSF) of 
the telescope. 

6.1.1. Lens mass profiles 

We use elliptically-symmetric distributions with 
power-law profiles to model the dimensionless surface 
mass density of the lens galaxies, 

where 7' is the radial power-law slope (with 7' = 2 
corresponding to isothermal), ^e is the Einstein radius, 
and q is the axis ratio of the elliptical isodensity con- 
tours. Various studies have shown that the power- 
law profile provides accurate descriptions of lens galax- 
ies (e.g.. [Cavazzj et al , 2007; H umphrev fc Buotc 20l3; 
Koopmans et al Ilmi Kuger. et a l. 2010; Barnabe; et al.l 
2011h . In particular. iSuvu et al, (20091 ) found that 
the grid-based lens potential corrections from power-law 
models were only ~ 2% for B1608-I-656 with interact- 
ing lens galaxies, thus validating the use of the sim- 
ple power-law models even for complicated lenses. We 
note that the surface brightness of the main deflector 
in RXJ1131— 1231 shows no signs of interaction (Section 
I6.1.2P and it is therefore much simpler than the case of 
B1608-f 656, further justifying the use of a simple power- 
law model to describe the mass distribution within the 
multiple images. 



Cosmological constraints from time-delay lenses 



7 



The Einstein radius in Equation (|12[) corresponds to 
the geometric radius of the critical curveQ and the mass 
enclosed within the isodensity contour with the geometric 
Einstein radius is 



(13) 



that depends only on a robust quantity in lensing. 

The deflection angle and lens potential of the power- 
law profile are computed following Barkana (1998). For 
each lens galaxy, the distribution is suitably translated 
to the position of the lens galaxy and rotated by the 
position angle cf) of the lens galaxy (where (/> is a free 
parameter, measured counterclockwise from 62). Since 
the satellite galaxy is small in extent, we approximate 
its mass distribution as a spherical isothermal mass dis- 
tribution with 7g = 2 and = 1 in Equation (jl2p . The 
(very small) impact of the satellite on cosmographic in- 
ferences is discussed in Section [5751 

Our coordinate system is defined such that 61 and 62 
point to the west and north, respectively. The origin of 
the coordinates is at the bottom-left corner of the ACS 
image containing 160x160 pixels. 

In addition to the lens galaxies, we include a constant 
external shear of the following form in polar coordinates 
6 and ip: 



V'oxt(^', </5) = -70x16*^ COS 2(</J- (/)oxt), 



(14) 



where 7ext is the shear strength and i/icxt is the shear 
angle. The shear position angle of (/)oxt = 0° corresponds 
to a shearing along the 61 direction whereas (pcxt — 90° 
corresponds to a shearing in the 62 direction. 

We do not include the external convergence Kext at this 
stage, since this parameter is completely degenerate with 
-Da* in the ACS and time-delay modeling. Rather, we 
use 1)^°'^''' = (1 — Kext)-DAt for the lensing and time-delay 
data, and information on Koxt will come from kinematics 
and lens environment in Section [7] to allow us to infer 
Dai- 

6.1.2. Lens light 

For the light distribution of the lens galaxies, we use 
elliptical Sersic profiles, 



I {61,92) = Aexp 



•cff 



(15) 

where A is the amplitude, k is a constant such that i?cff 
is the effective (half-light) radius, (7l is the axis ratio, and 
'T'sorsic is the Sersic index (jScrsic 196^). The distribution 
is suitably rotated by positions angle and translated 
to the galaxy positions (0i,l, ^2,l)- We find that a single 
Sersic profile for the primary le ns galaxy leads to signi f- 
icant residuals, as was found bv lClaeskens et al.l ()20Q6D . 
Instead, we use two Sersic profiles with common cen- 
troids and position angles to describe the lens galaxy G. 
For the small satellite galaxy that illuminates only a few 

The critical curve of Kpi in Equation II12I I is symmetric about 
di and 02 J and the geometric radius is -^Skm^^^hi^, where 6iong 
(^short) is the distance of the furthest (closest) point on the critical 
curve from the origin. 



pixels, we use a circular Sersic profile with rigcrsic = 1- 
This simplifying assumption has no effect on the mass 
modeling since the light of the satellite is central and 
compact, and thus does not affect the light or mass of 
the other components. 

6.1.3. Source light 

To describe the surfac e brig h tness distribution of the 
lensed source, we follow ISuvui (f20ll and use a hybrid 
model comprised of (1) point images for the lensed AGNs 
on the image plane, and (2) a regular grid of source sur- 
face brightness pixels for the spatially extended AGN 
host galaxy. Modeling the AGN point images inde- 
pendently accommodates variations in the fluxes arising 
from microlensing, time delays and substructures. Each 
AGN image therefore has three parameters: a position in 
61 and O2 and an amplitude. We collectively denote these 
AGN parameters as v. The ext e nded source on a grid 
is modeled following ISuvu et"aLl (j2Q06f ). with curvature 
regularization. 

6.1.4. PSF 

A PSF is needed to model the light of the lens galax- 
ies and the lensed source. We use stars in the field to 
approximate the PSF, which has been shown to work 
suffic iently well in modehng galaxy-scale strong lense s 
(e.g.. [Marshall et al.l[2007l: ISuvu et aLll^OOl ISuvull20Tl) . 
In particular, we adopt the star that is located at 2.'4 
northwest of the lens system as the model of the PSF. 

6.1.5. Image pixel uncertainties 

The comprehensive mass and light model described 
above captures the large-scale features of the data very 
well. However, small-scale features in the image might 
cause misfits which, if not taken into account, may lead 
to an underestimation of par amete r unce rtainties and 
biased parameter estimates. ISuvul (|2012D found that 
by boosting the pixel uncertainty of the image surface 
brightness, the lens model parameters can be faithfully 
recovered with realistic estimation of uncertainties. 

Following this study, we introduce two terms to de- 
scribe the variance of the intensity at pixel i of the ACS 
image ^acs, 



2 _ 2 

pix,i bkgd 



/C?ACS/ 



(16) 



where crbkgd is the background uncertainty, / is a scal- 
ing factor, and dp^cs,i is the image intensity. The second 
term, /dACS.ij corresponds to a scaled version of Poisson 
noise for the astrophysical sources. We measure txbkgd 
from a blank region in the image without astrophysical 
sources. We set the value of / such that the reduced 
is ~1 for the lensed image reconstruction (see, e.g., 
ISuvu et al] 120061 for details on the computation of the 
reduced that takes into account the regularization on 
the source pixels). Equation ([T5)) by design downweights 
regions of high intensities where the residuals are typ- 
ically most prominent. This allows the lens model to 
fit to the overall structure of the data instead of reduc- 
ing high residuals at a few locations at the expense of 
poorer fits to the large-scale lensing features. The resid- 
uals near the AGN image positions are particularly high 
due to the high intensities and slight saturations in some 
of the images. Therefore, we set the uncertainty on the 
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inner pixels of the AGN images to a very large number 
that effectively leads to these pixels being discarded. We 
discard only a small region in fitting the AGN light, and 
increase the region to minimize AGN residuals when us- 
ing the extended source features to constrain the lens 
mass parameters. 

6.2. Likelihoods 

The model-predicted image pixel surface brightness 
can be written as a vector 



"■ACS 



Bg + BLs 



Nag 

E 



N 



(17) 



where B is a blurring operator that accounts for the PSF 
convolution, g is a vector of image pixel intensities of the 
Sersic profiles for the lens galaxy light, L is the lensing 
operator that maps source intensity to the image plane 
based on the deflection angles computed from the param- 
eters of the lens mass distributions (such as 7', 6'e, 7cxt), 
s is a vector of source-plane pixel intensities, A'agn(= 4) 
is the number of AGN images, and ai(i>') is the vector of 
image pixel intensities for PSF-convolved image i of the 
AGN. 

The likelihood of the ACS data with image pixels 

is 



-P(dACs|7'7^'E,7oxt,i7) 
= y"dsP(dAcs|7'.^E,7oxt,rj,s)P(s), 



(18) 



where 



-P(dACs|7',^'E,7cxt,T7, s) 



JVagn 

n 



i=l 

1 



ACS.i 



27rcr, 



exp 



2<^pix,, 

\9, - e 



P|2 



2<jf 



(19) 



In the first term of this likelihood function, is just the 
normalization 



Zd = (27r 



piX,2 ; 



(20) 



f^ACS.i is the surface brightness of pixel i, 
(i^(-.g j(7', 7cxt, ''7, •s) is the corresponding predicted 
value given by Equation (|17p (recall that t] arc the 
remaining lens model parameters to which the ACS data 
are sensitive) , and a^^^ ^ is the pixel uncertainty given 
by Equation (|16|) . The second term in the likelihood 
accounts for the positions of the AGN images, modeled 
as independent points (i.e. non-pixelated sources) in 
the image. In this term, 6i is the measured image 
position (listed in Table [1]), ai is the estimated posi- 
tional uncertainty of 0.005", and 0f (7', 7cxt, ''7) is 
the predicted image position given the lens parameters. 
(Notice that this second term does not contribute to 
the marginalization integral of Equation p^ .) The 
form of P{s) for the source intensity pixels and the 
resulting analytic expression for the marginali zation in 
Equation ([18]) are detailed in iSuvu et all (|2006f) . 



TABLE 1 
Lens model parameters 



Description 


Marginalized 
Parameter or optimized 
constraints 


Timc-dclay distance (Mpc) 


r)modcl 
^At 


1883: 


-85 


Lens mass distribution 


Centroid of G in 9i (arcsec) 


dl.G 


4.420 


+U.003 
-0.002 


Centroid of G in 82 (arcsec) 


d2,G 


3.932 


+0.004 
-0.003 


Axis ratio of G 


QG 


0.763 


+0.005 
-0.008 


Position angle of G (°) 


<t>G 


115.8 


+0.5 
-0.5 


Einstein radius of G (arcsec) 




1 64+°°^ 


Radial slope of G 


7' 


1 05+0.05 


Centroid of S in 81 (arcsec) 




4.323 




Centroid of S in 82 (arcsec) 


6*2,3 


4.546 




Einstein radius of S (arcsec) 


9e.s 




External shear strength 


7cxt 


0.089 


+0.006 
-0.006 


External shear angle (°) 


0cxt 


92tl 





Lens light as Sersic profiles 



Centroid of G in 81 (arcsec) ^l.GL 

Centroid of G in 82 (arcsec) S2,GL 

Position angle of G (°) </'GL 

Axis ratio of Gl 9GL1 

Amplitude of Gl ^GLi 

Effective radius Gl (arcsec) ^?cfI,GLi 

Index of Gl risoisicGLi 

Axis ratio of G2 <?GL2 

Amplitude of G2 ^GL2 

Effective radius of G2 (arcsec) ficff,GL2 
Index of G2 

^sersic. GL2 

Centroid of S in 81 (arcsec) 6*1, SL 

Centroid of S in 82 (arcsec) S2,SL 

Axis ratio of S gsL 

Amplitude of S Agj^ 

Effective radius of S (arcsec) R^g gL 
Index of S 



4 .,,+U.OUi 
4 +0.001 

121.6±;^j 
n 070+0.004 

O.»'»_o.o03 

2 49+0 01 
^■^^-0.01 

n qs+0 03 
'j-yj_o.o3 
0.84910:004 

0.89to.03 

S62+0 009 

U.OOZ_Q QQg 

1 p;q+003 
^■"^^-0.03 
4.323 
4.546 

= 1 
34.11 
0.01 
1 



Lensed AGN light 


Position of image A 


in 


6\ (arcsec) 


dl.A 


2.383 


Position of image A 


in 


62 (arcsec) 


e2,A 


3.412 


Amplitude of image 


A 




O-A 


1466 


Position of image B 


in 


di (arcsec) 


^l.B 


2.344 


Position of image B 


in 


82 (arcsec) 


^2,B 


4.594 


Amplitude of image 


B 




as 


1220 


Position of image C 


in 


61 (arcsec) 


9l,C 


2.960 


Position of image C 


in 


62 (arcsec) 


d2,C 


2.300 


Amplitude of image 


C 




a-C 


502 


Position of image D 


in 


81 (arcsec) 


^l,D 


5.494 


Position of image D 


in 


82 (arcsec) 


^2,D 


4.288 


Amplitude of image 


D 






129 



Notes. There are a total of 39 parameters that are optimized or 
sampled. Parameters that are optimized are held fixed in the sam- 
pling of the full parameter space and have no uncertainties tabu- 
lated. Changes in these optimized parameters have little effect on 
the key parameters for cosmology (such as 1)^°'^°'). The tabulated 
values for the sampled parameters are the marginalized constraints 
with uncertainties given by the 16**^ and 84*'^ percentiles (to indi- 
cate the 68% credible interval). For the lens light, two Sersic pro- 
files with common centroid and position angle arc used to describe 
the primary lens galaxy G. They are denoted as Gl and G2 above. 
The position angles are measured counterclockwise from positive 
62 (north). The source surface brightness of the AGN host is mod- 
eled on a grid of pixels; these pixel parameters (s) are analytically 
marginalized and are thus not listed. 
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The likelihood for the time delays is given by 



n 



1 



exp 



(21) 



where Ati is the measured time delay with uncer- 
tainty (JAt.i for image pair i=AB, CB, or DB, and 
Atf (i:)^^°'^°\7',6'E,7cxt,»7) is the corresponding pre- 
dicted time delay computed via Equation ^ given the 
lens mass model parameters. 

The joint likelihood for the ACS and time 
delay data that appears in Equation (|T0)) . 
PidACS,At\DAt,l',OE,Jc^t,'n,i^e^t), is just thc 
product of the likelihoods in Equations and ([?!]) . 
We assign uniform priors over reasonable ranges for all 

mo(' 
At 



the lens parameters: _D^°^''', 7', ^e, 7cxt and ij. 



6.3. MCMC sampling 

We model the ACS image and time delays with 
GleeH- The ACS image has 160x160 surface bright- 
ness pixels as constraints. There are a total of 39 lens 
model parameters that are summarized in Table [TJ This 
is the most comprehensive lens model of RXJ1131— 1231 
to date. The list of parameters excludes the source sur- 
face brightness pixel parameters, s, which are analyti- 
cally marginalized in computing the likelihood (see, e.g., 
ISuvu fc Halkolall2010l for details) . With such a large pa- 
rameter space, we sequentially sample individual parts of 
the parameter space first to get good starting positions 
near the peak of the PDF before sampling the full pa- 
rameter space. The aim is to obtain a robust PDF for 
the key lens parameters for cosmography: 1)^°^''', 7', 6'e, 
and 7ext- 

For an initial model of the lens light, we create an an- 
nular mask for the lensed arc and use the image pixels 
outside the annular mask to optimize the lens Sersic pro- 
files. The parameters for the light of the satellite is fixed 
to these optimized values for the remainder of thc analy- 
sis since the satellite light has negligible effect on 
and other lens parameters. Furthermore, we fix the ccn- 
troid of the satellite's mass distribution to its observed 
light centroid. We obtain an initial mass model for the 
lenses using the image positions of the multiple knots in 
the so urce that are identified following IBrewer &: LewisI 
(j2008[) . Specifically, we optimize for the parameters that 
minimize the separation between the identified image po- 
sitions and the predicted image positions from the mass 
model. We then optimize the AGN light together with 
the light of the extended source while keeping the lens 
light and lens mass model fixed. The AGN light parame- 
ters are then held fixed to these optimized values. Having 
obtained initial values for all the lens model parameters 
to describe the ACS data, we then proceed to sample the 
lens parameters listed in Table [T] using a MCMC method. 
In particular, we simultaneously vary the following pa- 
rameters: all mass parameters of G, the Einstein radius 
of S, external shear, the extended source intensity dis- 
tribution, the lens light profile of G. Glee employs sev- 

a len s modeling software developed by S. H. Suyu and 
A. Halkola IISuvu fc Halkolall2010t) based on lSuvu et"all (20m) and 
IHalkola et al.l (12003) 



eral of the methods of iDunklev et al.l (120051 ) for efficient 
MCMC sampling and for assessing chain convergence. 

6.4. Constraints on the lens model parameters 

We explore various parameter values for the AGN light 
and the satellite Sersic light, try different PSF mod- 
els, and consider different masks for the lensed arcs and 
the AGNs (which are fixed in the MCMC sampling). 
These variations have negligible effect on the sampling 
of the other lens parameters. The only attribute that 
changes the PDF of the parameters significantly is the 
number of source pixels, or equivalently, the source pixel 
size. We try a series of source resolution from coarse 
to fine, and the parameter constraints stabilize start- 
ing at ~50x50 source pixels, corresponding to source 
pixel sizes of ~0'.'05. Nonetheless, the parameter con- 
straints for different source resolutions are shifted sig- 
nificantly from one another. Different source pixeliza- 
tions minimize the image residuals in different manners, 
and predict different relative thickness of the arcs that 
provides in formation on the lens profile slope 7' (e.g., 
ISuvul[20T^ . To quantify this systematic uncertainty, we 
consider the following set of source resolutions: 50x50, 
52x52, 54x54, 56x56, 58x58. 60x60, and 64x64. The 
likelihood P(dACS, At\D^f^\ 7', 0E, 7cxt), which is pro- 
portional to thc marginalized posterior of these parame- 
ters P(£'At"^°''7''^E,7oxt|dACS, At) since thc priors are 
uniform, is plotted in Figure [3] for each of these source 
resolutions. The scatter in constraints among the vari- 
ous source resolutions allows us to quantify the system- 
atic uncertainty. In particular, we weight each choice of 
the source resolution equally, and combine the Markov 
chains together. In Table [TJ we list the marginalized pa- 
rameters from the combined samples. 

We show the most probable image and source recon- 
struction for the 64x64 resolution in Figure Only the 
image intensity pixels within the annular mask shown in 
thc top-middle panel are used to reconstruct the source 
that is shown in the bottom-right panel. A comparison 
of the top-left and bottom-left panels shows that our lens 
model reproduces the global features of the ACS image. 
The time delays are also reproduced by the model: for 
the various source resolutions, the is ~2 for the three 
delays relative to image B. There are some small resid- 
ual features in the bottom- middle panel of Figure |4l and 
these cause the shifts in the parameter constraints seen 
in Figure [3] for different source pixel sizes. The recon- 
structed host galaxy of the AGN in the bottom-right 
panel shows a compact central peak, which is proba- 
bly thc bulge of the spiral source galaxy, embedded in 
a more diffuse patch of light (the disk) with knots/spiral 
features. 

6.5. Understanding the external shear 

The inferred external shear is 7cxt = 0.089 ± 0.006 
(marginalizing over all other model parameters) from 
modeling the ACS image and the time delays. The ex- 
ternal shear may provide information on the amount of 
external convergence, since they originate from the same 
external structures. However, the high 7ext found in our 
model could potentially be attributed to deviations of 
thc primary lens from its elliptical power-law descrip- 
tion; if this were the case, some of 7ext would in fact 
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Fig. 3. — Posterior of the key lens model parameters for the lensing and time-delay data. Each color represents a particular source 
resolution that is the dominant source of systematic uncertainty in the modeling of the ACS data. The solid curves are a Gaussian fit to 
the PDF by weighting each source resolution chain equally. The contours/shades mark the 68.3%, 95.4%, 99.7% credible regions. 



be internal shear. To gauge whether the modeled shear 
is truly external, we also considered a model that in- 
cludes a constant external convergence gradient. This 
introduces two additional parameters: k' (gradient) and 
4>K, (the position angle of the gradient, where 0^ = 
corresponds to positive k gradient along the positive 62 
direction, i.e., north). The ACS data allow us to con- 
strain k' = (S.ltg:;^) X IQ-^arcsec^i and 0^ = 87 ± 2°. 
The convergence gradient is aligned along the same di- 
rection as the external shear within 5° and has a sensible 
magnitude, suggesting that the shear is in fact truly ex- 
ternal, and is likely due to mass structures to the east of 
the lens. 

To investigate the origin of the external shear, we 
construct a wide-field R-band image from the COSMO- 
GRAIL monitoring images that is shown in Figure [S) 
The lens system is indicated by the box, and the galax- 
ies (stars) in the field are marked by solid (dashed) 
circle s, identified using SExtractor (jBertin fc ArnoutsI 
[19961) . Overlaid on the i mage within th e dash ed box 
are the X-ray contours from lChartas et ahl (j2009() . show- 
ing the presence of a galaxy cluster at z = 0.1 that 
is located at 158" northe ast of the lens () Morgan et alj 
120061: iChartas et aI]|2009D. Using the measured X-ray 
luminositv in lChartas et all (|2009f ). we estimate that the 
contribution of the cluster to the external shear at the 



lens is only a few percent. Nonetheless, large-scale struc- 
tures associated with the cluster and the plethora of mass 
structures to the east of the lens could generate addi- 
tional shear. The fact that our modeled external shear 
and convergence gradients both point toward mass struc- 
tures in the east that are visible in Figure [5] is a further 
indication that the modeled shear is indeed external. We 
will use this external shear in Section [7. 2 1 to constrain the 
external convergence. 

6.6. Propagating the lens model forward 

To facilitate the sampling and marginalization of the 
posterior of the cosmological parameters in Equations 
pUj) and pT|) . we approximate the overall likelihood of 
c'acs and Ai from the multiple source resolutions in Fig- 
ure [3] with a multivariate Gaussian distribution for the 
interesting parameters 7',6'E,7oxt and marginal- 
izing over the nuisance parameters rj. This approxima- 
tion allows the value of P(<iAcs, Ai|7', 6'e, 7oxt, £'a°'^''') 
to be computed at any position in this 4-dimensional pa- 
rameter space. Note that in contrast to the other param- 
eters, the Einstein radius of the primary, 9e, is well deter- 
mined, with minimal degeneracy with other parameters. 
This robust quantity is used in the dynamics modeling of 
the lens galaxy. The approximated Gaussian likelihood 
provides an easy way to combine with the stellar kine- 
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Fig. 4. — ACS image reconstruction of the most probable model with a source grid of 64x64 pixels. Top left: observed ACS F814W 
image. Top middle: predicted Icnscd image of the background AGN host galaxy. Top right: predicted light of the lensed AGNs and the lens 
galaxies. Bottom left: predicted image from all components, which is a sum of the top-middle and top-right panels. Bottom middle: image 
residual, normalized by the estimated la uncertainty of each pixel. Bottom right: the reconstructed source. Our lens model reproduces 
the global features of the data. 



matics and lens environment information for measuring 

DAt- 

7. CONSTRAINING Koxt 

In this section, we fold in additional information on the 
lens galaxy stellar kinematics and density environment to 
constrain the nuisance parameter Kcxt- 

7.1. Stellar kinematics 

We follow ISuYU et all (|201CI[ ) and model the velocity 
dispersion of the stars in the primary lens galaxy G, 
highlighting the main steps. The three-dimensional mass 
density distribution of the lens galaxy can be expressed 
as 



Pair) = (Kext - 



Til) 



1 



rr 



(22) 

Note that the projected mass of the lens galaxy en- 
closed within is (1 — Kext)AfE, while the projected 
mass associated with the external convergence is Kgxt-^E; 
the sum of the two is the Einstein mass Me that was 
given in Equation p3| . We employ spherical Jean's 
modeling to infer the line-of-sight velocity dispersion, 
g^fTT, 7^ 0E, rani, KpTtt) i pom pg by assuming the Hern- 



quist profile (Hernquistlll99Cll) for the stellar distr ibution 
(e.s.. IBinnev fc TremaindfToSTt ISuvu elall [201(11) 1^ An 

ISuvu et al.l 1I20T0I ') found that IHernguistl llT990h and Ija^ 
1)19831 ) stellar distribution functions led to nearly identical cosmo- 



anisotropy radius of rani = corresponds to pure radial 
stellar orbits, while rani — > 00 corresponds to isotropic 
orbits with equal radial and tangential velocity disper- 
sions. We note that cr^ is independent of Hq, but is 
dependent on the other cosmological parameters (e.g.,w 
and f2do) through Scrit and the physical scale radius of 
the stellar distribution. 

The likelihood for the velocity dispersion is 



P((T|7r,7',6'E,rani,Koxt) 



exp 



(cr - cr^(7r,7',6'E,ra„i,K;ext)) 



21 



2a2 



(23) 



where a = 323kms~^ and Co- = 20kms~^ from Sec- 
tion l4.3l Recall that the priors on 7' and 9e were assigned 
to be uniform in the lens modeling. We also impose a 
uniform prior on rani in the range of [0.5, 5]i?off for the 
kinematics modeling, where the effective radius based on 
the two-component Sersic profiles in Table[T]is 1'.'85 from 
the photometrvF^ The uncertainty in i?cff has negligible 
impact on the predicted velocity dispersion. The prior 
PDF for TT is discussed in Section ISTTl while the PDF for 
/vext is described in the next section. 

7.2. Lens environment 



logical constraints. 

Before unblinding, we used an effective radius of 3'.' 2 based 
on a single Sersic fit. The larger R^ff changes the inference of D/\t 
at the < 0.5% level. 
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Fig. 5. — 11.5' X 10.5' R-band image obtained from stacking 60 hours of the best-quaUty images in the COSMOGRAIL monitoring. The 
lens system is marked by the box near the center. Galaxies (stars) in the field are indicate d by solid (dashed) circles. The radius of the 
solid circle is proportional to the flux of the galaxy. X-ray map from IChartas et al.l l|2009l l are overlaid on the image within the dashed 
box. The concentrations of mass structures to the cast of the lens are consistent with the modeled external shear and convergence gradient 
directions. 

We combine the relative galaxy counts from Sec- 
tion |4]4l the measured external s hear in Sectionl6.4[ a nd 
the Millennium Simulation (MS; ISpringel et al.ll2005D to 



obtain an estimate of P{Ko-xt ^Rnv . 7Rxt . MS) . T his builds 
on the approach presented in lSuvu et al.l (|2010( ) that used 
only the relative galaxy counts. 

Tracing rays thr ough the Millennium Simulation (see 
iHilbert et aLll2009l for details of the method), we create 
64 simulated survey fields, each of sohd angle 4x4 deg^. 
In each field we map the convergence and shear to 
the source redshift Zg, and catalog the gala xy content, 
which we derive from the galaxy model by iGuo et all 
()2010D . For each line of sight in each simulated field, 
we record the convergence, shear, and relative galaxy 
counts in a 45" aperture having I-band magnitudes be- 
tween 18.5 and 24.5. These provide samples for the PDF 
Kext: 7ext J c^env|MS). Wc assume that the constructed 
PDF is app l icable to strong- lens lines of sight, following 
iSuvu et all (|201Clf ) who showed that the distribution of 



Kcxt for strong lens lines of sight are very similar to that 
for all lines of sight. 

Structures in front of the lens distort the time delays 
and the images of the lens/source, while structures be- 
hind the lens further affect the time delays and images of 
the source. However, to model simultaneously the mass 
distributions of the strong lens galaxies and all struc- 
tures along the line of sight is well beyond current ca- 
pabilities. In practice, the modeling of the strong lens 
galaxies is performed separately from the description of 
line-of-sight structures, and we approximate the effects 
of the lines-of-sight structures into the single correction 
term Kc^t, whose statistical properties we estimate from 
the Millennium Simulation. 

By selecting the lines of sight in the Millennium Sim- 
ulation that match the properties of RXJ1131 — 1231, we 
can obtain -P(denv|'«oxt, 7oxt, MS) /'(^oxt) and simultane- 
ously marginalize over 7ext in Equation pop . We as- 
sumed a uniform prior for j^xt in the lensing analysis. 
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Fig. 6. — The effective prior probability distribution for the ex- 
ternal convergence Koxt from combining ray tracing through the 
Millennium Simulation with the relative galaxy counts of 1.4 and 
the modeled external shear of 0.089 it 0.006. Dashed line: the 
convergence distribution for all lines of sight; Dot-dashed line: the 
convergence distribution for lines of sight with relative galaxy count 
rir = 1.4 it 0.05; Dotted line: the convergence distribution for all 
lines of sight weighted by the likelihood for 7oxt from the lens 
model; Solid line: the 7cxt-weightcd convergence distribution for 
lines of sight with rtr = 1.4 it 0.05. The effective prior for Kcxt used 
in the final cosmological parameter inference is described by this, 
most informative, distribution. 

such that P(7oxt) is a constant. The lensing likelihood is 
the only other term that depends on 7oxt, and from Sec- 
tion l6.4[ the lensing likelihood provides a tight constraint 
on 7ext that is approximately Gaussian: 0.089 ± 0.006. 
We can therefore simplify part of Equation ([TO)) to 

/ d7ext-P('iACS, Atl^At, 7') ^E, 7oxt, Kext) 
•P(de„v|Kext, 7ext, MS) 
- P{d ACS, ^t\D At, l',dE,Keyit) 

•P(denv|Kcxt, 7cxt = 0.089 ± 0.006, MS), (24) 

where the above approximation, i.e., neglecting the co- 
variance between 7ext and the other parameters in the 
lensing likelihood and then marginalizing 7oxt separately, 
is conservative since we would gain in precision by includ- 
ing the covariances with other parameters. Furthermore, 
by Bayes' rule, 

P(de„v|Kcxt, 7oxt - 0.089 ± 0.006, MS)P(Atext) 

«P(Kcxt|de„v,7cxt = 0.089 ±0.006, MS), (25) 

which is precisely the PDF of n^xt by selecting the sam- 
ples in P(Koxt, 7oxt, c^envlMS) that satisfies dgnv with a 
relative galaxy count within 1.4±0.05, and subsequently 
weighting these samples by the Gaussian likelihood for 
7cxt- This effective prior PDF for Koxt that is con- 
structed from the weighted samples, P(Koxt|'ienv, 7cxt = 
0.089 ± 0.006, MS), is shown by the solid line in Figure 

m 

8. TIME-DELAY DISTANCE OF RXJ1131-1231 

We combine all the PDFs obtained in the previous sec- 
tions to infer the time-delay distance Da*- 

8.1. Cosmological priors 

As written above, we could infer the time delay dis- 
tance Z^At directly, given a uniform prior. However, we 
are primarily interested in the cosmological information 
contained in such a distance measurement, so prefer to 



infer these directly. The posterior probability distribu- 
tion on Z?At can then be obtained by first calculating the 
posterior PDF of the cosmological parameters tt through 
the marginalizations in Equations (jlip and (jlOp . and 
then changing variables to P'a*- Such transformations 
are of course straightforward when working with sam- 
pled PDFs. 

In Table [2j we consider the following five cosmological 
world models, each with its own prior PDF P(7r): 

• UHq: Uniform prior PDF for between and 
150 kms~^Mpc~^ in a ACDM cosmology with 
51a = 1 — ~ 0.73. This is similar to the typical 
priors that were assumed in most of the early lens- 
ing studies, which sought to constrain Hq at fixed 
cosmology. 

• UuiCDM: Uniform priors on the parameters 
{iJo,fidc7w} in a fiat wCDM cosmology, where w 
is time- independent and flm = 1 — ^^dc- 

• WMAPTwCDM: The prior PDF for the parameters 
{Hq, J7do7 w} is taken to be the post erior PDF from 
the WMAP 7-year data set (Koma tsu et al.ll201lD . 
assuming a fiat wCDM cosmology, where w is time- 
independent and fini = 1 — ^^do- 

• WMAP7oACDM: The prior PDF for the parame- 
ters {Hq, VI a, r^k} is taken to be the posterior PDF 
from the WMAP 7-year data set, assuming an open 
(or rather, non-flat) cosmology, with dark energy 
described by A and ilk = 1 ^ ^^A — as the cur- 
vature parameter. 

• WMAP7owCDM: The prior PDF for the parame- 
ters {Hq, fide, w, ilk} is taken to be the posterior 
PDF from the WMAP 7-year data set, assum- 
ing an open wGDM cosmology, where w is time- 
independent and $lk = l^^^do~^^m is the curvature 
parameter. 

8.2. Posterior sampling 

We sample the posterior PDF by weighting samples 
drawn from the prior PDF with t he joint li kelihood func- 
tion evaluated at those points (jSuvu et al . 2010). We 
generate samples of the cosmological parameters tt from 
the priors listed in Table [2] We then join these to sam- 
ples of Kcxt drawn from P(Kcxt) from Section 17.21 and 
shown in Figure [6l and to uniformly distributed samples 
of 7' within [1.5,2.5] and rani within [0.5,5] Poff- Rather 
than generating samples of from the uniform prior, we 
obtain samples of 9e directly from the Gaussian approx- 
imation to the lensing and time-delay likelihood since 
is quite independent of other model parameters (as 
shown in Figure [3]). This boosts sampling efficiency, and 
the 9^ samples are only used to evaluate the kinematics 
likelihood. 

For each sample of {tt, Koxt, 7', ^ani, ^e}, we obtain the 
weight (or importance) as follows: (1) we determine P'a* 
from 77 via Equation ([2]), (2) we calculate Daj"*"' via 
Equation 0, (3) we evaluate P{dAcs, ^t\D'^f''\y) 
based on the Gaussian approximation shown in Figure [3] 
for Dlf'=^ and 7', (4) we compute P(c7|7r, 7', /text, 6'e, 
^ani) via Equation (j23p . and (5) we weight the sam- 
ple by the product of P(dAcs, At|I?S°'i°i, 7') and 
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TABLE 2 

Priors on cosmological parameters 



Prior 



Description 



UwCDM 



uniform Hq in [0, 150] km s ^ Mpc ^ , 

= 0.27, Ua = 0.73, w = -1 
uniform Hq in [0, 150] kms~^ Mpc~^, 
uniform = 1 ~ in [0, 1], 
uniform w in [—2.5,0.5] 
WMAP7i«CDM WMAP7 for {Ho, U^^, w} with flatness 

and time-independent w 
WMAP7oACDM WMAP7 for {Ha, Ha, Un^} with w = -1 
in open ACDM 

WMAP7o«)CDM WMAP7 for {Hq, ^de, w, Q^} with time- 
independent w and curvature Qy^ 

P(cr|7r, 7', Kcxt, ^E, ''ani) from the previous two steps. 
The projection of these weighted samples onto tt or D^t 
effectively marginalizes over the other parameters. 

8.3. Blind analysis in action 

As a brief illustration of our blind analysis approach, 
we show in the left panel of Figure [7] the blinded plot of 
the time-delay distance measurement. For all cosmolog- 
ical parameters such as Dai, £'a°'^°'; ^0, w, ftni, etc., 
we always plotted their probability distribution with re- 
spect to the median during the blind analysis. There- 
fore, we could use the shape of the PDFs to check our 
analysis without introducing experimenter bias by blind- 
ing the absolute parameter values. When we marginal- 
ized the parameters during the blind phase, our analysis 
code also returned the constraints with respect to the 
median. For example, the blinded time-delay distance 
for the WMAPTwCDJVI cosmology would be O^Hq Mpc. 
We used this particular cosmology as our fiducial world 
model during the blind analysis. In the remainder of the 
paper, we show the unblinded results of RXJ1131— 1231. 
The comparison with the gravitational lens B1608-I-656 
and other cosmological probes was performed after we 
unblinded the analysis of RXJ1131— 1231; otherwise, the 
blind analysis would be spoiled by such a comparison 
since the results of these previous studies were already 
known. 

8.4. Posterior PDF for Dai 

We show in the right-hand panel of Figure [7] the un- 
blinded probability distribution of the time-delay dis- 
tance for the first three cosmological models in Table O 
The priors, shown in dotted lines, are broad and rather 
uninformative. When including information from dACS, 
At, cr, and Koxt, we obtain posterior PDFs of DAt 
for RXJ1131— 1231 that are nearly independent of the 
prior, demonstrating that time-delay lenses provide ro- 
bust measurements of DAt- We find that the data con- 
strain the DAt to RXJ1131— 1231 with ~6% precision. 

We can compress these results by approximating the 
posterior PDF for DAt as a shifted log normal distribu- 
tion: 



P{DAt\Hn,nd,,w,n„ 
1 

■ exp 



(log(a; - Ap) - md)^ 
2al 



,(26) 



V27r(a; - AD)crD 
where x = £'At/(lMpc), Ad = 1425.6, /zd = 6.4993 



and ctd = 0.19377. This approximation accurately re- 
produces the cosmological inference, in that Hq is recov- 
ered within < 1% in terms of its median, 16"^ and 84"^ 
percentile values for the WMAP7 cosmologies we have 
considered. The robust constraint on Da* serves as the 
basis for cosmological inferences in Section HI 

8.5. Sources of uncertainty 

Our DAt measurement accounts for all known sources 
of uncertainty that we have summarized in Table [3] The 
dominant sources are the first three items. The precision 
for the time delay is the la uncertainty as a fraction of 
the measured value for the longest delay, AtoB- For the 
lens mass model and line-of-sight contributions, we de- 
fine the precision as half the difference between the 16"^ 
and 84"^ percentiles of the PDF for D'^"'^"^ from Sec- 
tion [6]4] in fractions of its median value and for Kext from 
Section [7.2l in fractions of 1, respectively. The remaining 
sources of uncertainty are collectively denoted by "other 
sources" , and the two main contributors to this category 
are the peculiar velocity of the lens and the impact of 
the satellite. 

Spectroscopic studies of the field of RXJ1131— 1231 
indicate that the lens is in a galaxy group with a ve- 
locity dispersion of 429193^ km s'^ (|Wong et all 1201 If ). 
RXJ1131— 1231 is the brightest red sequence galaxy in 
this group, and is thus likely to be near the center of 
mass of the group halo with a small peculiar vcloc- 
ity relative to the group ([Za bludoff & Mulchacy 199^; 
IWilliams et all [2001 IGeorge et al.ll2012l) . However, the 
group could be moving relative to the Hubble flow due to 
nearby large-scale structures. The one-dimensional rms 
galaxy pecuhar velocity is typically < 300 km s~^ (e.g., 
iMosauera fc Kochanek(l201lt lPeebieslll993[ ). A peculiar 
velocity of 300kms~^ for RXJ1131 — 1231 would cause 
DAt to change by 0.8%I3 A similar peculiar velocity 
for the lensed source has a much smaller impact on DAt , 
changing it by only 0.2%. We note that the peculiar 
velocities of lenses are stochastic, and this source of un- 
certainty should average out in a sample of lenses. 

We have explicitly included the satellite in our lens 
mass model in Section [51 However, there is some de- 
generacy in apportioning the mass between the satellite 
and the primary lens galaxy since lensing is mostly sensi- 
tive to the total mass enclosed within the lensing critical 
curves (approximately traced by the arcs). The more 
massive the satellite, the less massive the primary lens 
galaxy. Owing to its central location and the degener- 
acy with the mass of the main deflector, we expect the 
impact of the satellite on the difference in gravitational 
potential between the multiple images to be very small. 

To assess the effect of the mass of the satellite on our 
DAt inference, we consider an extreme model where the 
satellite has zero mass. In this case, we require a more 
massive primary lens galaxy with higher 9^ to fit the 
lensing features, as expected. The resulting 1?^°'^'^' and 
y from this model are consistent with that of the orig- 
inal model, but with larger parameter uncertainties due 
to poorer fits without the satellite. Even if we use the 
overestimated of the primary lens from this extreme 

T he change in red shift due to pecuhar velocities is described 
in, e.e:., IHarrisonl l l 19741 ). 
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Fig. 7. — Blinded (left) and unblinded (right) PDFs for D^t, showing the RXJ1131 — 1231 posterior constraints on Da* (solid) given 
assorted priors for the cosmological parameters (dotted, labeled). See the text and Table [2] for a full description of the various priors. 
RXJ1131— 1231 provides tight constraints on D/^t, which translates into information about Qm, f^dei a^nd particularly Hq. 



TABLE 3 

Error budget on time-delay distance of RX,J1131 — 1231 



Description 


uncertainty 


time delays 


1.6% 


lens mass model 


4.6% 


line-of-sight contribution 


4.6% 


other sources 


<1% 


Total (Gaussian approximation) 


6.7% 


Total (full sampling) 


6.0% 



Notes. The other sources of uncertainty that contribute at the 
<1% include the peculiar veloci ty of the lens and the impact of the 
satellite. Details are in Section 18.51 The Gaussian approximation 
simply adds the uncertainties in quadrature, providing a crude es- 
timate for the total uncertainty based on the full sampling of the 
non-Gaussian PDFs. 

model for the kinematics, we find that the effect on the 
inferred Dai is at the < 1% level. 

In Table 131 we list the total uncertainty of 6.7% based 
on a simple Gaussian approximation where we add up the 
uncertainties of each contribution in quadrature. This is 
close to the more accurate 6.0% based on proper sam- 
pling that takes into account the non-Gaussian distribu- 
tion (e.g., of Koxt) and the inclusion of the stellar ve- 
locity dispersion. We note that the sampling does not 
include explicitly the other sources that contribute at 
the <1% level; however, they are practically insignif- 
icant in the overall error budget. Most of the uncer- 
tainty in Z?At comes from the lens mass model and the 
line-of-sight contribution. To reduce the uncertainty on 
RXJ1131— 1231's Da* would require a better model of 
the source intensity distribution that depends less sensi- 
tively on the source pixel size (po ssibly via an adaptive 
sourc e pixelization scheme (e.g., IVegetti fc KoopmansI 
l2009f )). and a better characterization of Koxt by using 
more observational information from the field. Investiga- 
tions are in progress to improve Kcxt constraints (Greene 
et al. in preparation, CoUett et al. in preparation). 



9. COSMOLOGICAL INFERENCE 

We now present our inference on the parameters of the 
expanding Universe and compare our results to other cos- 
mographic probes. Specifically, our D^t measurement 
for RXJ1131— 1231 provides information on cosmology 
that is illustrated in Section 19.11 We compare the re- 
sults to that of B1608+656 to check for consistency in 
Section before combining the two lenses together in 
Section [931 We then compare the constraints from the 
two time-delay lenses to a few other cosmological probes 
in Section 19.41 

9.1. Constraints from RXJ1131- 1231 

We have seen that the RXJ1131 — 1231 D^t measure- 
ment is nearly independent of assumptions about the 
background cosmology. While Da* is primarily sensitive 
to iJo, information from Da* must be shared with other 
cosmological parameters via the combination of angular 
diameter distances. Therefore, cosmological parameter 
constraints will depend somewhat on our assumptions 
for the background cosmology. In this section we con- 
sider the first three cosmologies listed in Tabled Ui?o, 
UwCDM, and WMAPTwCDM. 

With all other parameters fixed in the U7?o cosmol- 
ogy except for Hq, all our knowledge of Da* is con- 
verted to information on Hq. We therefore obtain a 
precise measurement of Hq = 78.7il^4'g kms~^ Mpc~^ for 
RXJ1131-1231 with a 5.5% uncertainty 

Next, we relax our assumptions on r^do, ^^m and w, 
and consider the UwCDM and WMAPTwCDM cosmolo- 
gies in a flat Universe. Figure |51 shows the result- 
ing constraints. The contours for the UwCDM cos- 
mology with vertical bands in the panels illustrates 
that the time-delay distance is mostly sensitive to Hq. 
The constraint on Hq breaks the parameter degeneracies 
in the WMAP7 data set, and we obtain the following 
joint parameter constraints for RXJ1131— 1231 in com- 
bination with WMAP7: Hq = 80.0l^;^ kms'^ Mpc-\ 
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Fig. 8.— RXJ1131-1231 marginalized posterior PDF for Ho, 
fide and to in flat uiCDM cosmological models. Contours/shades 
mark the 68.3%, 95.4%, 99.7% credible regions. The three sets of 
contours/shades correspond to three different prior /data set com- 
binations. Shaded red: RXJ1131 — 1231 constraints given by the 
UuiCDM prior; dashed blue: the prior provided by the WMAP 
seven-year data set alone; solid black: the joint constraints from 
combining WMAP and RXJ1131-1231. 



do 



0.79 ±0.03, and w 



1 9c;+0.17 



9.2. Comparison between RXJl 131- 1231 and 
B 1608+ 656 

How do the results of RXJ1131 — 1231 compare with 
that of B1608+656? We show in Figure [9] the overlay 
of the cosmological constraints of RXJ1131 — 1231 and 
B1608+656 in VHq (top panel) and UiuCDM (bottom 
panel). To investigate the consistency of the two data 
sets, we need to consider their likelihood functions in 
the multi-dimensional cosmological parameter space: in- 
consistency is defined by ins ufficient overlap betwe en the 
two likehhoods. We follow [Marshall et al.l ()2006f ). and 
compute the Bayes Factor F in favor of a single set of 
cosmological parameters and a simultaneous fit: 



RtB\ 



F : 



(27) 



are the likelihoods of the 



where and 
RXJ1131-1231 and B1608+656 data respectively, com- 
puted at each prior sample point. See the Appendix for 
the derivation of this result. 

For the cosmology VHq, the Bayes Factor is 3.2; for 
UuiCDM, it takes the value 3.8. For comparison, with 
two one-dimensional Gaussian PDFs, F takes the value 
of 1 when the two distributions overlap at their 2a 
points, and is about 3.6 when they overlap at their la 
point. From this we conclude that the results from 
RXJ1131— 1231 and B1608+656 are consistent with each 
other. We do not detect any significant residual sys- 
tematics given the current uncertainties in our measure- 
ments. 
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Fig. 9.— Comparison of RXJ1131-1231 (solid red) with 
B1608-I-656 (dotted blue) in VHq (top) and UuiCDM (bottom) 
cosmologies. The two distributions overlap within 2cr. The cosmo- 
logical constraints from the two lenses are statistically consistent 
with each other: the ratio of the probability that the two lenses 
share global cosmological parameters to the probability that the 
lenses require independent cosmologies is 3.2 in UHq and 3.8 in 
UtoCDM. 

9.3. RXJ1131-1231 and B1608+656 in unison 

Having shown that RXJ1131-1231 and B1608+656 
are consistent with each other, we proceed to combine 
the results from these two lenses for cosmological in- 
ferences. In particular, we consider the constraints in 
the WMAP7wCDM and WMAP7oACDM cosmologies 
in Table [2 

We show in Figure[TU]the cosmological constraints from 
individual lenses in combination with WMAP7, and the 
combination of both lenses and WMAP 7. By combining 



Cosmological constraints from time-delay lenses 



17 



TABLE 4 

Cosmological constraints from RXJllSl — 1231 and 
B1608+656 IN combination with WMAP7 



WMAPTuCDM 



Cosmology 


Parameter 


Marginalized 


Precision 






Value (68% CI) 






Ho 


75.2+1- 


5.7% 


luCDM 






2.5% 




W 


_1 14+0-17 


18% 




Ho 


'■^-^-3.6 


4.0% 


oACDM 


Ua 


n 75+0-01 

U- 'O_o.02 


1.9% 




Ok 


U.UUCi_Q QQg 


0.6% 



Notes. The Hq values are in units of kms~i Mpc~i. The "preci- 
sion" in the fourth column is defined as half the 68% confidence 
interval, as a percentage of 75 for Hq, 1 for Q^^t ^Ai ^^nd f^ki a-nd 
— 1.0 for w. 

the two lenses, we tighten the constraints on Hq, f^dc and 
i7k- The precision on w does not improve appreciably. 
With its low lens redshift, RXJllSl — 1231 provide very 
little information to w in addition to that obtained from 
B1608-I-656. In Table |4l we summarize the constraints 
from the two lenses. 

9.4. Comparison of lenses and other cosmographic 
probes 

How do the robust time-delay distances from the strong 
lenses compare to the distance measures of other probes? 
We show in Figure [11] a comparison of the cosmologi- 
cal constraints of the two lenses, Baryon Acoustic Os- 
cillat i ons (BAO; e . g.. iPercival et al.l 120101: iBlake et al.) 
20111: iMeh ta et al.l 120 1^, and sup ernovae (SN; e.g. 



Hicken et a l. 2009: Suz uki et al.ll2012l ). when each is com- 
bined with WMAP7 in the owCDM cosmology. The 
figures are qualitative since the samples for WMAP7 
chain in the owCDM cosmology are sparse and we 
have smoothed the contours after importance sampling. 
Nonetheless, we see that the sizes of the contours are 
comparable, suggesting that even a small sample of time- 
delay lenses is a powerful probe of cosmology. Both the 
lenses and BAO are strong in constraining the curvature 
of the Universe, while SN provides more information on 
the dark energy equation of state. Lenses are thus highly 
complementary to other cosmographic p robes, partic- 
ularly the CMB an d SN (see also, e.g., iLindeil [20Tll: 
iDas fc LindeH[20T^ . Each probe is consistent with flat 
ACDM: = and w = -1 are within the 95% credible 
regions. 

In Figure 1121 we compare the precisions on ilk and 
w in oiuCDM from the following cosmological probes in 
com bination with WM AP: the Cepheids distance lad- 
der (jRiess et al.l |2011|F1. BAO from the S loan Digi- 
tal Sky Survey fSDS S) JPercival et al.l[2010l) . SN from 
the Union 2.1 sample (jSuzuki et al. 20121) . reconstr ucted 
BAO using the SDSS galaxies ([MehtaeLalJ ISQIJ) and 
our two time-delay lenses. We note that the precisions 
on the Cepheids and the time-delay lenses are only ap- 
proximates since the samples of WMAP 7 are sparse in 

i'^ To derive the constraints on fi^ a-nd ^ from the combination 
of Cepheids and WMAP7, we importance sample the WMAP7 
chain by a Gaussian likelihood centered on Hq = 73.8 — 1.475(io + 
1) km s~l Mpc~i with a width of 2.4kms~i Mpc~i IIRiess et al.l 
120111) . The -1.475(io + 1) correspond s to t he tilt in the Hq-w 
plane shown in Figure 10 of lHiess et al.l 120111) . 
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Fig. 10. — Cosmological constraints from the combination of 
RXJ1131-1231 and B1608-I-656 assuming WMAP7toCDM (top) 
and WMAP7oACDM (bottom) cosmologies. The combined 
posterior PDF is shown by the solid contours, the PDF for 
RXJ1131-1231 in combination with WMAP7 is the shaded con- 
tours, and the PDF for B1608-I-656 in combination with WMAP7 
is the dotted contours. Contours/shades mark the 68.3%, 95.4%, 
99.7% credible regions. 
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Fig. 12. — Precision of cosmological constraints on VL^ and 
w for five probes each in combination with WMAP7 in an 



otuCDM cosmology: SDSS BAO UPercival et al.l |20TO| 'I. the two 
time-delay lenses RXJ1131-1231 and B1608-j-656 (this work), 
SN ISiSu ki ct al. 20 lJ), Cepheid s l IRiess et al.|[2MlD . and recon- 
structed BAO HMehta et al.ll2012l ). Precision for and w is de- 
fined as half the 68% CI as a percentage of 1. 



owCDM due to the large parameter space. Nonethe- 
less, the histogram plot shows that time-delay lenses are 
a valuable probe, especially in constraining the spatial 
curvature of the Universe. 

10. SUMMARY 

We have performed a blind analysis of the time-delay 
lens RXJ1131 — 1231, modeling its high precision time de- 
lays from the COSMOGRAIL collaboration, deep EST 
imaging, newly measured lens velocity dispersion, and 
mass contribution from line-of-sight structures. The data 
sets were combined probabilistically in a joint analysis, 
via a comprehensive model of the lens system consisting 
of the light of the source AGN and its host galaxy, the 
light and mass of the lens galaxies, and structures along 
the line of sight characterized by external convergence 
and shear parameters. The resulting time-delay distance 
measurement for the lens allows us to infer cosmologi- 
cal constraints. From this study, we draw the following 
conclusions: 



1. Our comprehensive lens model reproduces the 
global features of the HST image and the time 
delays. We quantify the uncertainty due to the 
deflector gravitational potential on the time-delay 
distance to be at the 4.6% level. 

2. Based on the external shear strength from the lens 
model and the overdensity of galaxy count around 
the lens, we obtained a PDF for the external con- 
vergence by ray tracing through the Millennium 
Simulation. This Koxt PDF contributes to the un- 
certainty on Z?At also at the 4.6% level. 

3. Our robust time-delay distance measurement of 6% 
takes into account all sources of known statistical 
and systematic uncertainty. We provide a fitting 
formula to describe the PDF of the time-delay dis- 
tance that can be used to combine with any other 
independent cosmological probe. 

4. The time-delay distance of RXJ1131-1231 is 
mostly sensitive to _ffo, especially given the low red- 
shift of the lens. 

5. Assuming a flat ACDM with fixed 17a = 0.73 
and uniform prior on Hq. our unblinded 
i/o measurement from RXJ1131— 1231 is 

Mpc"^ 



78.7t^;^kms-i^ 



6. The constraint on helps break parameter 
degeneracies in the CMB data. In combina- 
tion with WMAP7 in wCDM, we find Eq = 
80.0t^;?kms-iMpc"\ 17dc = 0.79 ± 0.03, and 

w = — 1.251q2i. These are statistically consis- 
tent with the results from the gravitational lens 
B1608+656. There are no significant residual sys- 
tematics detected in our method based on this com- 
bined analysis of the two systems. 



7. By combining RXJ1131-1231, B1608-h656 and 
WMAP7, we derive the following constraints: 
ffo = 75.21^2 km s-iMpc-\ r^dc = ^■'^^'^-nm 
and w = -1.14;q2o in flat wCDM, and Eq = 



73.ltlikm.s- 
^■^O^toml in open ACDM. 



iMpc-i, riA ^ 0.75tom and n 
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8. Our measurement of the Hubble constant is com- 
pletely independent of those b ased on the local 
distance ladder method (e.g., iRiess et all 120111 : 
iFreedman et al.l [2012), providing an important 
consistency check of the standard cosmological 
model and of general relativity. 

9. A comparison of the lenses and other cosmological 
probes that are each combined with WMAP7 shows 
that the constraints from the lenses are compara- 
ble in precision to various state-of-the-art probes. 
Lenses are particularly powerful in measuring the 
spatial curvature of the universe, and are comple- 
mentary to other cosmological probes. 

Thanks to t he dedicated monitoring by the COSMO- 
GRAIL (e.g.. IVuissoz etldl [2008t ICourbin et all [20lTI: 
iTewes et alll2012bllal ) and lKochanek et al.l (|2006r collab- 
orations, the number of lenses with accurate and pre- 
cise time delays are increasing. Deep HST imaging for 
three of these lenses will be obtained in cycle 20 to allow 
accurate lens mass modeling that turns the delays into 
distances. Current and upcoming telescopes and sur- 
veys including the Panoramic Survey Telescope & Rapid 
Response System, Hyper-Suprime Camera on the Sub- 
aru Telescope, and Dark Energy Survey expect to detect 
hundreds of AGN lenses with dozens of delays measured 
(|Oguri fc Marshalll 1201(1 . Ultimately, the Large Syn- 
optic Survey Telescope will discover thousands of time- 
delay lenses, painting a bright future for cosmography 
with gravitational lens time delays. 
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APPENDIX 

QUANTIFYING DATA SET CONSISTENCY VIA THE BAYES FACTOR B 

iMarshall et al.l (j2006f ) invite us to consider the following two hypotheses: (1) Hs'°'^''', in which the two lenses share 
a common set of cosmological parameters tt — {Ho, fido, w}, and (2) H'"'^, in which each of the two lenses is provided 
with its own independent set of cosmological parameters, tt^ = {H^, il^^, w^} and tt^ = {Hq, i^f,,, w^}, with which 
to fit the data. Each set of parameters covers the same prior volume as in Hs^°''^'. If the two data sets are highly 
inconsistent, only H'"*^ will provide a good fit to both data sets in a joint analysis. The question is, do the data require 
H''^'^, or is Hsi°'^^i sufficient? 

We quantify the answer to this question with the evidence ratio, or Bayes Factor, in favor of H^'"'''*': 



P(d",d^|Hg'°'^^') 



''' ^ . (Al) 



where we have collectively denoted all the data sets of RXJ1131— 1231 as and of B1608+656 as . Each of the 
terms on the right-hand side of the above equation can be written in terms of a multi-dimensional integral over the 
cosmological parameters. For example, we have (starting with the simpler terms in the denominators) 

P(d^|H''^'i) = /d37r'^P(d«|7r^,H'"^)P(7r^|H'"'^), (A2) 



where P(d^|7r^ , H''^'*) is the likelihood of the RXJ1131-1231 data sets (the weights for the cosmological samples) that 
we denote by L^ . Equation (|A2p is then just the ensemble average of the samples' likelihood values, 

P(£Z^|H'"'*) = (L^). (A3) 
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For P(d^|H'"<^), we can rewrite the lik elihood P(d^|7r^, H'"'^) in terms of to make use of P{D^^\d^, H'""!) given 
by Equation (35) of iSuvu et ahl (|20Tol) : 

Pid\rr,H ) = PiDl^) • ^^^^ 

The ratio = P(d^|H'"'i)/P(i:>|j|H'"'i) is a constant factor since the prior on D^^. is uniform; thus, we obtain 

P(d^|H'"'*) = Z^(L^), (A5) 

where L is given by the likelihood of the time-delay distance P{Df^{ir^)\d^ , H'"'*). 
Finally, for the numerator in Equation (|A1[) . we have 

P{d^, d^|Hsi°ba^) = J d^Ti- P(d^|7r, H^'"''''') P(d^|7T, H«'°''''') P(7r|Hs'°''^') 

= Z^{L^L^), (A6) 

where the constant is the same as that in P(<i^|H'"'^) since the parameterization of the cosmology for each inde- 
pendent lens is identical to that of the global cosmology (i.e., tt^ and tt are the same cosmological parameterization). 
Substituting Equations (|A3|) . (jASp . and (|A6p into Equation ()A1|) . we obtain 

(iR)(pB)' ^"^'^ 
which can be readily computed given the values of and (the weights) that we have for each cosmological sample. 



